Melting of graphene: from two to one dimension 
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Abstract 

The high temperature behaviour of graphene is studied by atomistic simulations based on an 
accurate interatomic potential for carbon. We find that clustering of Stone- Wales defects and 
formation of octagons are the first steps in the process of melting which proceeds via the formation 
of carbon chains. The molten state forms a three-dimensional network of entangled chains rather 
than a simple liquid. The melting temperature estimated from the two-dimensional Lindemann 
criterion and from extrapolation of our simulation for different heating rates is about 4900 K. 
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The structural properties of graphene are as exceptional as its electronic structure 1 -. Be- 
ing the simplest two-dimensional (2D) membrane, graphene is a natural benchmark of our 
understanding of physics in 2D^. Melting in 2D is usually described in terms of creation of 
topological defects, like unbound disclinations and dislocations 4 . In the hexagonal lattice of 
graphene, typical disclinations are pentagons (5) and heptagons (7) while dislocations are 
5-7 pairs. In carbon systems, pairs of 5-7 dislocations that may be created by rotation of a 
carbon bond are called Stone- Wales (SW) defects. SW defects are non-topological defects 
since they have Burger vector and Frank vector equal to zerc£ and thus do not introduce 
any torsion or curvature in the system. Although graphene is a 2D crystal, it bears also 
many analogies to other graphitic structures such as fullerenes and nanotubes for which the 
melting has been studied in Refs||] and 78 respectively. Since it is known that melting of 
small particles is essentially different from that of bulk crystals, the relation between the 
melting mechanisms in fullerenes and nanotubes and that in graphene is not a priory known. 
Besides, fullerenes possess an intrinsic curvature that might drastically affect the whole pic- 
ture of melting^. Our atomistic simulations give a scenario of the melting of graphene as 
the decomposition of the 2D crystal in a 3D network of ID chains. The molten phase is 
similar to the one found for fullerenes^. Similar structures can be seen in simulations of 
nanotubes^. We find that, for graphene, a crucial role is played by SW defects. It is the 
clustering of SW defects that triggers the formation of octagons which are the precursors 
for the spontaneous melting around T m « 4900 K in our simulations. This temperature 
is higher than the value of 4000 K found for fullerenes- and close to the 4800 K found for 
nanotubes^, but not as high as the value 5800 K estimated for graphene by extrapolation to 
infinite radius nanotubes in Ref s). 

We study the melting of graphene by Monte Carlo simulations in the NVT and NPT 
ensembles (constant number of particles N, constant volume V or constant pressure P and 
constant temperature T) with periodic boundary conditions in the plane for samples of 
N = 1008, iV = 4032 and N = 16128 atoms. The interatomic interactions are calculated 
with the LCBOPII potential that we have shown to describe well the elastic and thermal 
properties of graphene 1 ^ as well as the liquid phase 1 ^. It is important to note that this 
potential allows bond breaking and formation with realistic energy barriers. Also it gives a 



much better description of the lattice dynamics - than the Tersoff potential used in Ref. 



7. To 



speed up equilibration and improve the sampling of long wavelength out-of-plane distortions, 
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we also used collective wave moves as described in Ref . [13|. More details on the computational 
procedure can be found in our previous works^ii^. To estimate the melting temperature 
T m we have performed simulations where the temperature is increased at slower and slower 
heating rates as we will describe. 

In Fig. [TJwe show a typical configuration of graphene on the way to melting at 5000 K > 
T m . The coexistence of crystalline and melted regions indicates a first order phase transition. 
The most noticeable features are the puddles of graphene that have melted into chains. The 
melted areas are surrounded by disordered 5-7 clusters, resulting from the clustering and 
distortion of SW defects. They have the smallest formation energy and start appearing at 
~ 3800 K. In Fig. [TJ also isolated and pairs of SW defects are present whereas we never 
observe isolated pentagons, heptagons or 5-7 dislocations. Possibly, the less energetically 
favourable ring structures observed in the simulations of nanotubes^ are due to a less accurate 
description of the energetics given by the Tersoff potential. In comparison with the previous 
studies for fullerenes and nanotubes, the much larger size of our samples allows to investigate 
the role of the spontaneous formation and annihilation of SW defects in the melting process. 
In Refill the effect of SW defects on the melting temperature was studied by adding SW 
defects by hand. 

In Fig. |2]we show the average number of heptagons^ as a function of temperature together 
with the Arrhenius behaviour with the formation energy E$w — 4.6 eV given by LCBOPII^. 
We see that only for T < 4000 K the fit is quite good which demonstrates that all heptagons 
are part of almost isolated SW defects as also confirmed by visual inspection of typical 
configurations. At higher temperatures, deviations indicating a higher formation energy 
per heptagon signals the observed clustering of SW. We find that 5-7 clusters act as nuclei 
for the melting. This is in contrast with graphite where melting is initiated by interplanar 
covalent bond formation. Close inspection shows that regions with 5-7 clusters favor the 
transformation of three hexagons into two pentagons and one octagon (Fig. [3]) that we never 
see occurring in the regular hexagonal lattice far from the 5-7 clusters or near isolated SW 
defects. The octagons, in turn, are the precursors for the formation of larger rings. Due 
to the weakening of the bonds forming the relatively small angles in the pentagons around 
them, the atoms forming these larger rings tend to detach from the lattice and form chains. 

When melting is completed the carbon chains form an entangled 3D network with a 
substantial amount of three-fold coordinated atoms, linking the chains. Therefore, this low 
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density structure reminds rather a polymer gel than a simple liquid. In fact the radial 
distribution function shown in the inset of Fig. [1] displays a very sharp peak at distances 
smaller than the interatomic distance in graphene, due to the formation of shorter double 
bonds in the chains. In typical simple liquids, instead, the first peak shifts to larger values of 
interatomic distances when going from a crystal to a liquid 15 . This high temperature phase 
has also been found for molten fullerenes^. 

In Fig. H] we show the evolution of the potential energy and structural properties of 
an initially flat, graphene layer as a function of Monte Carlo steps (1 step is equal to iV 
displacement trials) at T = 4750 K (below melting) and at T = 5000 K where melting 
occurs within about 10 7 Monte Carlo steps around step 2.5 x 10 7 . We see that the rise of the 
potential energy at melting, is mirrored by the growth of the number of chains (nc) at the 
expenses of the six-member rings (R6) of the crystalline phase. The number of eight-member 
rings (R8) that are formed close to clusters of SW defects, instead, increases when melting 
starts and decreases when chains are formed, illustrating the melting mechanisms described 
previously. 

The Lindemann criterion^ 1 ^ v/ (u 2 ) m /d ~ 0.2 (where (u 2 ) m is the mean-square atomic 
displacement at the melting temperature and d is the interatomic distance at T = 0) is 
commonly used to estimate the melting temperature in 3D systems. Since in 2D the mean 
square displacement (u 2 ) is divergent at finite temperatures^ 7 -, we need to consider differences 



of atomic displacements^ 19 . Adapting the melting criterion used in Refjl8lto the honeycomb 
lattice of graphene we define the average quantity 



n 
i 



(1) 



where a = 1/y/Wpo where po is the 2D particle density at T = K, r« is the position of 
the i-th atom and where the sum over j runs over the n atoms closest to atom i. In Fig. [5] 
we show 7 n for n = 3, 9 and 12, namely, including one, two and three coordination spheres 
(see the inset in Fig. [5]). Since the difference between second and third neighbor distances 
is relatively small, the distinction between second and third neighbours becomes fuzzy at 
very high temperatures as one can see from the merging of the related peaks in the radial 
distribution in the inset of Fig. [TJ Therefore we think that 712 is more meaningful than 79. 
Interestingly, as shown in Fig. [5j melting occurs when 7 12 ~ 0.1 as found for the strictly 
2D triangular lattice in Ref.— . Instead, 73 remains much smaller due to the rigidity of the 



covalent nearest-neighbour bonds. 

To estimate the melting temperature we consider the effect of heating the sample at 
different rates. In Fig. we show the potential energy per atom (Fig. Eh) and 712 (Fig. Eb) 
for several heating rates, each half of the previous one. One can see that the temperature 
at which the energy and 712 suddenly grow, signaling the melting, moves from about 5200 
K towards the left saturating slightly above 4900 K. 

The closest system to graphene is graphite. The melting temperature of graphite has 
been extensively studied experimentally at pressures around 10 GPa and the results present 
a large spread between 4000 K and 5000 K 20 . With LCBOPII, free energy calculations give 
T m = 4250 K, almost independent of pressure between 1 and 20 GPaM: At zero pressure, 
however, graphite sublimates before melting at 3000 K— . Monte Carlo simulations with 
LCBOPII at zero pressure show that, at 3000 K, graphite sublimates through detachment 
of the graphene layers^. The melting of graphene in vacuum that we have studied here can 
be thought of as the last step in the thermal decomposition of graphite, the 2D graphene 
layers melting into a 3D liquid network of ID chains. Interestingly, formation of carbon 
chains has been observed in the melt zone of graphite under laser irradiation 22 Although 
the temperature T = 4900 K of spontaneous melting represents an upper limit for T m , our 
simulations suggest that T m of graphene at zero pressure is higher than that of graphite. 
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Figure 1: Top: Snapshot of graphene at T = 5000 K during the melting process in a NPT 
simulation [N = 1008, P = 0). All pentagons and heptagons are marked in red. Bottom: side 
view of the same sample when the melting is completed. Insert: radial distribution function (rdf) 
at T = 4750 K and T = 5000 K, below and above melting. After melting, the first peak is shifted 
to smaller distances, reflecting the chain formation and further structure, typical of the crystalline 
phase, is washed out. 
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Figure 2: Average number of heptagonal rings (R7) in the NVT ensemble as a function of temper- 
ature for the N = 16128 sample. The solid line shows the Arrhenius behaviour with the formation 
energy of SW defects given by LCBOPII, E$w = 4.6 eV which gives a good description only up to 
E ~ 4000 K whereas a fit to the data for T > 4400 (dashed line) gives a higher formation energy 
of tss 4.95 eV, due to SW clusterization. 



8 




9 




6000 



£ 4000 



2000 




3 50 



b) 



4750 K 
5000 K 



c) 

4750 K 

-5oooK m ;m 



1 



'I, 



500 

























-4750K 










5000 K 


* 

J! 1 , 


■ ■ ■ 1 





2 3 
MC steps 



x 10 



Figure 4: Behaviour of several quantities as a function of the number of Monte Carlo steps starting 
from a solid flat graphene layer (N = 16128) at T = 4750 K (below melting) and at T = 5000 
K where melting occurs within about 10 7 Monte Carlo steps around step 2.5 x 10 7 . From top to 
bottom: a) Total potential energy E in eV/atom; b) Number of six-member rings R6; c) Number 
of 8-member rings R8. Eight-member rings are formed close to clusters of SW defects. R8 increases 
when melting starts and decreases when chains are formed, illustrating the melting mechanism (see 
text); d) Number of chains nc with more than 3 connected two- fold coordinated atoms. 



10 




Q I I I I I I I 

1500 2000 2500 3000 3500 4000 4500 5000 

T, K 



Figure 5: Temperature dependence of 73,79 and 712 calculated including the closest three, nine 
and twelve neighbours (see inset). The last points indicated by arrows correspond to the onset of 
the liquid phase for which these quantities diverge for infinite systems. 
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Figure 6: Potential energy (a) and 712(b) as a function of temperature for graphene layer (N = 
1008) for five different heating rates, each two times slower than the previous one, namely 1, 1/2, 
1/4, 1/8, 1/16 degree Kelvin per 1000 Monte Carlo steps as indicated by the labels. 
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